Single and multi-trait genomic prediction for agronomic traits in Euterpe edulis

Popularly known as juçaizeiro, Euterpe edulis has been gaining prominence in the fruit growing sector and has demanded the development of superior genetic materials. Since it is a native species and still little studied, the application of more sophisticated techniques can result in higher gains with less time. Until now, there are no studies that apply genomic prediction for this crop, especially in multi-trait analysis. In this sense, this study aimed to apply new methods and breeding techniques for the juçaizeiro, to optimize this breeding program through the application of genomic prediction. This data consisted of 275 juçaizeiro genotypes from a population of Rio Novo do Sul-ES, Brazil. The genomic prediction was performed using the multi-trait (G-BLUP MT) and single-trait (G-BLUP ST) models and the selection of superior genotypes was based on a selection index. Similar results for predictive ability were observed for both models. However, the G-BLUP ST model provided greater selection gains when compared to the G-BLUP MT. For this reason, the genomic estimated breeding values (GEBVs) from the G-BLUP ST, were used to select the six superior genotypes (UFES.A.RN.390, UFES.A.RN.386, UFES.A.RN.080, UFES.A.RN.383, UFES.S.RN.098, and UFES.S.RN.093). This was intended to provide superior genetic materials for the development of seedlings and implantation of productive orchards, which will meet the demands of the productive, industrial and consumer market.


Introduction
Euterpe edulis, a palm tree native to the Atlantic Forest, is distributed along the entire length of this biome and in the gallery forests of the Brazilian cerrado [1][2][3]. It belongs to the genus Euterpe, which has become popular in recent years through açaí, processed pulp from the fruits of species of this genus [4]. The most exploited species is Euterpe oleracea, açaí, which was recognized as a fruit tree of economic expression in 2008 [5].
As a crop, Euterpe edulis has become popular under the name of juçaizeiro. It is a perennial species that can exceed 15 meters in height, with a single stem [6] that does not regrow [7]with extremely rare cases of tillers. Its inflorescences are infrafoliate, located in the transition a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Field experimental conduction
Experimental evaluations were carried out from 2018 to 2021 in the municipality of Rio Novo do Sul in the state of Espírito Santo, Brazil (Fig 1). The 275 genotypes evaluated were selected for good phytosanitary and physiological conditions and full reproductive age. The genotypes are located in a commercial plantation belonging to two companies that process juçaizeiro fruit pulp, Açaí juçara1 and Bonalotti1. The plantation was formed by a mixture of individuals that emerged by spontaneous development, through the action of dispersers, and the  by the owners by the sowing of seeds in the area. Commercial planting does not have pre-defined spacing and basically no productive management treatment, mowing is only done at harvest time. Consequently, the field evaluations followed a blank design, as the descent relationship between the genotypes was initially unknown.
Phenotyping of plants. Between 2018 and 2021, the number of bunches per plant (NB) was evaluated by visually counting the number of bunches with fruits of each matrix plant. In the years 2018, 2019 and 2021, the mass of fruits per bunch (MFB) (kg), rachis length (RL) (cm), equatorial diameter of fruits (EDF) (mm) and pulp yield (PY) (%) were evaluated. For the evaluation of MFB, in all evaluations, the same trained professional performed the determination of the harvest point of the bunches. The harvest of each genotype was carried out when the fruit maturation point reached the stage used for fruit processing in the industry. Each year, one bunch per plant was harvested, the fruits were separated from the rachillas, and they were weighed on a scale with an accuracy of 0.1 g. For RL, after separating the fruits, the length of the rachis of the inflorescence was measured with a tape measure.
A sample of fruits was taken from each genotype, packed in properly identified plastic bags and transported to the Plant Biometry laboratory at the Federal University of Espírito Santo, where morphometric evaluations of fruits and seeds were carried out in a completely randomized design. For EDF, data were measured in millimeters (mm) obtained using a 6" digital caliper (Zaasprecision1), performed on five fruits individually, as recommended by Marçal et al. (2015) [30], who observed high repeatability values for these characteristics. The authors concluded that five measurements are necessary to carry out the measurement to reach coefficients of determination of 95%. Measurements above this amount increase costs and evaluative time, bringing little additional information to the works.
The PY was estimated by the following relationship: In the equation, FFM is fruit fresh mass and SFM is seed fresh mass, measured by weighing four replicates of 25 fruits and seeds, using an analytical balance (0.0001g).
Genomic DNA extraction. Genomic DNA was obtained from leaf samples of the genotypes under study. Extraction was performed using the cetyltrimethylammonium bromide or CTAB method by Doyle (1990) [31] with modifications [32].
DNA concentrations and integrity were estimated using a NanodropTM 2000 spectrophotometer (Thermo Scientific). DNA quality was verified on 0.8% agarose gel. DNA genotypes prepared for genotyping using the DArTseqTM methodology were sent to the Service of Genetic Analysis for Agriculture (SAGA) in Mexico for high-throughput genotyping using the DArTseqTM technology The genome representation of the 275 genotypes was obtained from the reduction of DNA complexity using two restriction enzymes, HpaII and Msel. The ends of the cleaved fragments were linked to a code adapter and a common adapter to identify each sample. The fragments were amplified using PCR reaction; subsequently, equimolar amounts of amplification products from each sample of the 96-well microtiter plate were pooled, purified and quantified, then sequenced on the Illumina Novaseq 6000 System platform. The sequences were analyzed using Dartsoft14, an automated genomics data analysis program and DArTdb, a laboratory management system, developed and patented by DArT Pvt. Ltd. (Australia), generating SNP marker data as described by Killian et al. (2012) [33] and Sansaloni et al. (2020) [34].
Quality control of molecular markers. The dataset of codominant markers of the SNP type was submitted to quality control analysis in the R [35]. We removed the markers with Call Rate (CR) � 90% and Minor Allele Frequency (MAF) � 5%. After quality control, the marker dataset reduced by 81.75% from 44,457 markers to 8,112.
Phenotypic data analysis. The phenotypic values were corrected for the effect of years in the R software, according to the correction proposed by   [36]. The linear model used was: where y ik is the phenotypic value for the genotype i and the year k; μ is the population mean; a k is the fixed effect of the k th year; g i is the fixed effect of the i th genotype; and e ik is the random effect of residual, with e � Nð0; Is 2 e Þ. The empirical best linear unbiased estimates (eBLUEs) was calculated for each trait individually and their values were obtained byŷ i ¼ m þâ k þ g i .
Genomic prediction (GP). Genomic prediction was performed using the multi-trait (G-BLUP MT) and single-trait (G-BLUP ST) models. For this the SOMMER package version 3.4 [37] was implemented in the R software [35]. The G-BLUP MT model for the prediction of the genomic estimated breeding values (GEBVs) of the individuals used was: Where y is the vector of previously estimated eBLUES and structured as y ¼ ½ y 1 y 2 . . . y n � 0 ; y 1 ; y 2 . . . ; y n is the vector of observations for each characteristic; b is the vector of means of each characteristic structured as b ¼ ½ b 1 b 2 . . . b n � 0 and with X incidence matrix. u is vector of individual additive genomic genetic values of each trait structured as u ¼ ½ u 1 u 2 . . . u n � 0 . Z is the incidence matrix, with the variance structure given by u~N(0, G�S u ), where G is the genomic relationship matrix between individuals for additive effects, S u is the additive genetic covariance matrix and � denotes the Kronecker product; e is the random error vector with e~N(0, I�S e ) where S e is the residual covariance matrix. Covariance matrices can be written as: where s 2 u i and s 2 e i are, respectively, the additive and residual genetic variance associated with the i th trait with i = 1,. . .,n = 5; s 2 u ij and s 2 e ij are, respectively, the additive and residual genetic covariance associated with the i th and j th traits with j = 1,. . .,n = 5 and i 6 ¼ j. The variance and covariance components were obtained via the restricted maximum likelihood method (REML). The additive genomic relationship matrix (G) was obtained as described by VanRaden (2008) [38] by the centralization of the matrix of markers: The parameterization for the incidence matrix W is presented below and is in accordance with the classical theory of quantitative genetics [39]: where, p i and q i are allele frequencies of M and m, respectively.
With the genetic values of the individuals for the NB and MFB traits, the fruit production per plant (FPP) for each genotype was estimated.
The narrow-sense heritability of the i th trait was estimated following the equation below: whereŝ 2 u i andŝ 2 e i are the additive genetic variance and the residual variance of the i th trait, respectively. The genetic correlation between the i th and j th traits was obtained through the following equation: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffî where s 2 u i and s 2 u j are the additive genetic variance of the i th and j th trait, respectively. s 2 u ij is the additive genetic covariance between the i th and j th traits.
Predictive ability and cross-validation. The predictive ability (r y;û ) and the standard error of r y;û , were estimated through the cross-validation procedure, randomly subdividing the population into 5 folds. Thus, 220 genotypes were used for training set and 55 genotypes were used for validation set. For each fold, the r y;û was obtained by the correlation between the predicted GEBV's (û) and the corrected phenotypic values (y).
Genetic selection based on selection index. With the GEBVs predicted by the G-BLUP ST method, we initially estimated the production of fruit production per plant (PFP) by multiplying the GEBVs between NB and MFB. The selection of the best genotypes was based on multiple variables. For this, the method developed by Mulamba and Mock (1978) [40] was used, which is based on the sum of the individual ranks of each characteristic, creating a global rank. For the selection, the parameters of NB, MFB, PFP, RL and PY in the positive direction and EDF in the negative direction were considered. Two processes were performed before rank summing. The first was to transform the EDF values into classes, in order to create four classes. Class I for small fruits (values below the first quartile of the distribution of the genetic values); class II for small/medium size fruits (genotypes with values between the first and second quartile of the distribution of the genetic values); class III for medium fruits (genotypes with values between the second and third quartile of the genetic values), and class IV for large fruits (genotypes with values above the third quartile of the genetic values). The second change was the normalization of the number of ranks for each trait, to avoid traits with fewer ranks having a greater influence on the selection process. The standardization followed the expression below: where r ij is the standardized rank value for genotype i and characteristic j; nr j is the number of ranks of trait j and p ij and is the rank of genotype i for the trait j.
In order to compare the selective efficiency between the G-BLUP ST and G-BLUP MT models, it was estimated the expected selection gain for the population of the next selection cycle and for the commercial seed donor population. For both methods and all characteristics, with the exception of PFP, estimates were obtained by: where, SG is the expected selection gain for characteristic i by model j; h 2 ij is the heritability of characteristic i estimated by model j and SD i is the selection differential for characteristic i in model j estimated based on corrected phenotypic values.
The SG for PFP was estimated based on the product of the SG estimates of NB and MFB. Comparison between methodologies. Cohen's Kappa coefficient [41] was used to analyze the agreement between the best selected individuals between the G-BLUP ST e G-BLUP MT models, for the commercial seed donor population (formed by the six best ranked individuals), and the next selective cycle population (formed by the 50 best ranked individuals). Cohen's Kappa coefficient is given by: where NAO are the number of observed agreements, NAEC is the number of expected agreements by chance, and NOA is the number of analyzed observations [42].

Data description
A summary of the descriptive statistics including mean, standard deviation, and maximum and minimum values for the five agronomic traits evaluated in this work is shown in Table 1. It is noted that in general, the population presented a small variation in the average phenotypic response between the years evaluated.

Genetic parameters, genetic correlation and phenotypic correlation
The heritability and the residual and additive genetics variance components were estimated for G-BLUP ST and G-BLUP MT models (Table 2). Considering the standard error, all these estimates were statistically equivalent for all traits in the evaluated models. Heritabilities ranged from 0.29 to 0.80 for NB (GBLUP-MT) and PY (GBLUP-ST), respectively ( Table 2). The genetic and phenotypic correlations between the traits are shown in Fig 2, which were calculated using the (co)variance estimates obtained by the G-BLUP MT model, so that the genetic correlation estimates are presented on the upper diagonal and, on the lower diagonal, estimates of phenotypic correlations are presented.
Comparing the phenotypic and genotypic correlations between the traits in Fig 2, it is observed that only two correlation estimates had changes in their directions (RL and NB; MFB and NB), with the change in correlation between MFB and NB being more pronounced. The results obtained for the phenotypic and genotypic correlations can be classified from low to moderate [43]. For phenotypic correlation, the highest values observed were 0.41 (RL and MFB), in the positive direction, and -0.17 (RL and PY), in the negative direction. While, for the genotypic correlations between traits, the highest values were 0.40 (RL and EDF) in the positive direction, and -0.37 (NB and MFB) in the negative direction (Fig 2). E os erros padrões para as correlações genéticas, foram de magnitudes mais elevadas quando comparado aos estimados para as correlações fenotípicas. Standard errors for genetic correlations ranged from 0.18 (PY and EDF) to 0.37 (MFB and NB) and for phenotypic correlations, standard errors ranged from 0.01 (PY and EDF) to 0.12 (PY and RL) (S1 Table).

Predictive ability
The predictive ability (r y;û ) for the G-BLUP ST and G-BLUP MT models is shown in Table 3. For the G-BLUP ST, the r y;û ranged from 0.21 for NB to 0.49 for PY, and for the G-BLUP MT model, the r y;û ranged from 0.18 for NB to 0.48 for PY. In general, the values of r y;û observed for both methodologies were close, as well as NB and PY maintained their positions with the lowest and highest r y;û , respectively. The results show that both models are similar to perform the prediction of genomic genetic values.

Matrix selection and expected genetic advancement
The Cohen's Kappa coefficient was calculated for the top 20% of individuals (50) with the highest rank positions based on the Mulamba and Mock index [40], aimed to form the population of the next selection cycle. The agreement between the G-BLUP ST and G-BLUP MT models was 0.68 and is classified as substantial (0.60-0.80) [44] (Fig 3). However, when we reduce the selected population to six individuals, to build the commercial seed donor population, the coefficient is reduced to an agreement of 0.50, considered moderate. Even though the agreements of the Cohen's Kappa coefficient were classified from moderate (0.41-60) to substantial (0.61-80) [44], is observed in Fig 3, a great divergence in the ranking of individuals between the methods used. To support the comparison of the efficiency   between the methods, the selection gains (SG) estimate was obtained through the heritability and phenotypic means of the selected individuals, and the results are shown in Table 4. In general, higher SG resulting from the G-BLUP ST model is observed for both populations. The population mean, the phenotypic and genotypic behavior of the selected individuals, based on the information from the G-BLUP ST model and the field phenotypic, is shown in Fig 4. It is possible to observe that the phenotypic means of the six selected individuals, in the desired direction, have great performance in relation to the population average, being the improvement 31.13%, 51.05%, 98,40%, 11.56%, 15.57%, and 7.43% for NB, MFP, PFP, RL, PY and EDF, respectively.
The selection of the six best genotypes provides a change in phenotypic responses equivalent to 0.99 bunches, 1.70 kg, 10.43 kg, 7.08 cm, 5.09%, and -1.00 mm for the traits NB, MFP, PFP, RL, PY and EDF, respectively (Fig 4). With SG, it is expected an improvement on the averages in 0.41 bunches, -0.69 kg, 3.62 cm, 5.09%, and -0.47 mm, for NB, MFP, RL, PY, and EDF, respectively. For PFP, the genomic values were estimated by the product of the genomic values of MFB and NB, and the difference between the means of the selected and the general population was 3.71 kg (Fig 4).

Discussion
As there are few programs that are based on the improvement of juçaizeiro, there are also a few scientific works that are aimed at the application of selective techniques for individuals of this species. In this sense, the existing breeding programs are in the early stages of development, and thus, their techniques are basically based on classical methodologies with phenotypic information. In contrast to these processes, the present work innovates in being the first scientific study focused on genomic prediction of Euterpe edulis, with the objective of increasing selective accuracy and selection gains within the breeding cycle.

Genetic parameters, genetic correlation and phenotypic correlation
The h 2 u found in the present work for the traits under study indicate that the morphometric traits of fruits (EDF and PY) have a greater potential for heredity and genetic control, when compared to the other productive parameters related to clusters (RL, MFB and NB). This behavior was also observed in Euterpe oleracea (açaizeiro), a species of the same genus as juçaizeiro [45,46].
It is expected that traits with higher h 2 u , consequently, have a higher predictive ability [47]. The same was found by Legarra et al. (2008) [47], where the traits with the highest h 2 u in the Table 4

. Estimates of selection gains provided by the G-BLUP ST and G-BLUP MT models. Estimates of selection gains provided by the single-trait (G-BLUP ST) and
multi-trait (G-BLUP MT) models for the selected genotypes (S.I.) for the population of the next selection cycle (50 genotypes), and the commercial seed donor population (six genotypes). present study (EDF and PY) also had the higher predictive ability (

PLOS ONE
is known that for quantitative traits, several environmental factors can affect phenotypic behavior. In this sense, it can be assumed that the genetic control of NB, MFB and RL is reduced due to the exposure of these traits for a longer period of time to environmental effects and more variable conditions that are not controlled, for example, direction of insertion of the bunch into the matrix plant, exposure to inclement weather, feeding the fauna, among others. The development of fruits in the bunch will still occur, for example, in a fraction of the time of the infructescence development cycle, which will be located in the same region of the plant, with a similar influence on luminosity, temperature, and availability of nutrients. Therefore, genetic control can be expected to be greater.
Published scientific studies evaluating associations between Euterpe edulis traits are scarce in the literature [12,30,48]. However, these analyses are fundamental for understanding the interactive behavior between traits, revealing linear cause and effect responses that are essential for the agricultural development of the crop and for breeding programs, as it allows the determination of different types of practical strategies.
The knowledge of correlations in the genetic scope is important in breeding programs because they are heritable and can be used in indirect selection [49]. However, the evaluation of phenotypic behavior is also essential, since selection carried out indirectly based on genotypic information, and can lead to unwanted results in the field. This occurs when the phenotypic associations have the opposite direction to the genetic correlations considered for use in indirect selection.
The (co)variances estimated by the G-BLUP MT, made it possible to estimate the Pearson correlation for genetic and phenotypic effects between the traits RL, EDF, MFB, NB and PY (Fig 2), allowing a counterpoint between these estimates. In this sense, we noticed that practically all the associations preserved their influence on the behavioral response, with the exception of the NB-RL and NB-MFB pairs, which had an inversion of the sense of their estimates. In this condition, indirect selection is compromised by the fact that selection by genotypic responses may not be expressed in the desired sense in the environment. With this, we can conclude that the analyses of phenotypic and genotypic associations are complementary and make it possible to determine practical actions, using indirect processes to enable increased gains in breeding programs.
Evaluating biometric characteristics of juçaizeiro fruits, Marçal et al. (2015) [30] reported a positive genetic association between EDF and fresh seed mass (0.75). Knowing that the juçaizeiro fruit is mostly composed of seed, and the pulp makes up a small portion of it-increasing EDF can reduce the percentage of pulp yield. This fact is confirmed by the association observed in Fig 2, which shows inversely proportional associations between EDF and PY (-0.17 phenotypic; -0.21 genetic), that is, the increase in fruit size results in a reduction in the percentage of pulp per fruit, harmful to industrial processing.
The results observed by Oliveira et al. (2015) [48] for the genetic correlation in Euterpe edulis, are close to those observed in the present work between the characteristics NB, PY and RL. In both situations, the magnitudes of the associations can be classified from very weak to weak [43]. Oliveira et al. (2015) [48] reported associations of 0.20, -0.10 and -0.02 for NB-PY, PY-RL and NB-RL, respectively, while the present work obtained estimates of 0.07, -0.12 and -0.02, respectively. Corroborating these results, Farias Neto et al. (2016) [50] also found a weak correlation between NB-RL (-0.01) for Euterpe oleracea.
Regarding the correlation results observed in Fig 2, we can determine that the use of indirect selection among the evaluated characteristics would be inefficient. Therefore, the effect that the characteristics presented among themselves were mostly low, making the use of this practice unfeasible due to the small gains that could be obtained.

Predictive ability
In order to increase selection efficiency, choosing the best methodology for analyzing and predicting the genetic values of individuals is one of the fundamental aspects to be taken into account in a breeding program. In this study, we compare G-BLUP ST e G-BLUP MT models for the evaluation of genomic prediction in a base population of juçaizeiro. The existence of genetic correlations between the selection traits is the basis for the advantages presented by G-BLUP MT [51] and the absence of correlation could lead to equivalent or even superior results when using G-BLUP ST [51]. Therefore, we can conclude that, depending on the analyzed traits and the existence of correlations, it is expected that the use of the GBLUP MT model will increase the accuracy of the evaluation predictions [52].
However, the superiority of G-BLUP MT over G-BLUP ST for r y;û was not observed in the present work. The results obtained by the G-BLUP MT model showed similar values of r y;û when compared to the G-BLUP ST model, which may have been caused by the linear associations that varied from low to moderate intensity (� 0.42) (Fig 2) and, high standard of error associated with these estimates (S1 Table). As reported by Calus and Veerkamp (2011) [53], who observed an improvement in r y;û when marker information was added to the model for the traits that had associations less than 0.50. In this sense, it is expected that the flow of genetic and phenotypic information provides higher r y;û when compared to G-BLUP ST.
The slightly lower performance of G-BLUP MT for predictive abilities can be explained by the low genetic and residual correlations between traits, as reported by Runcie and Cheng (2019) [54], who attribute such results to the fact that low correlations can cause imprecision in the estimates of genetic and residual covariance parameters may result in reduced model performance. However, as the r y;û were close, there is no evidence of such negative effects on the covariance estimates, deducing that the flow of information between the characteristics was not sufficient to improve the r y;û of the model G-BLUP MT in relation to the G-BLUP ST, considering that marker information was used by both models.
The r y;û for GBLUP-ST and GBLUP-MT performed statistically equivalent for most traits, with the exception of the RL, which showed r y;û superior to G-BLUP-ST (0.29) compared to G-BLUP MT (0.22). The results diverge from patterns observed by other studies [27,51,52,[55][56][57][58][59], in which superiority is observed for the G-BLUP MT model. As it was designed to benefit from the existence of genetic correlations between traits, it is expected that the G-BLUP MT model will present better results than the G-BLUP ST [11,53]. Due to the similarity of the results to the predictive ability of the models G-BLUP ST and G-BLUP MT, it is not possible, just for this parameter, to indicate the best model for the selection of juçaizeiro genotypes.

Genotype selection and expected genetic advancement
To build the next selective cycle breeding population, and in order to maintain genetic variability, a selection intensity of approximately 20% was defined. In this condition, the Kappa agreement observed among the approaches evaluated was classified as substantial (0.68), considering that its value was between 0.60 and 0.80 [44]. This similarity between the approaches can be explained by the fact that as the variables under study present low to moderate correlations [43], both methods ended up leading to similar results.
We should add that the high Kappa agreement between the methods is also associated with the conditions of the phenotypic data set used in the study, which has a low degree of imbalance. Thus, the predictive abilities of the approaches did not have a great impact on the selection of different genotypes between models. However, the influence on the ranking of the best individuals was strongly impacted, as shown in Fig 3. When considering the six best genotypes, the selection of the best genotypes had greater divergence between the models, in which the Kappa value was 0.50, evidencing the need to determine the best model to carry out the selection process, aiming at greater future gains. In this way, the models are chosen according to their ability to predict more accurately the GEBVs of the analyzed individuals. However, in the present study, r y;û was not sufficient to determine the best method, due to similar estimates observed.
Due to the similarities of the r y;û and the great divergence in the ranking for the Kappa agreement, the choice of the best model for the selection of the juçaizeiro genotypes was based on the SG. In order to standardize the comparison between the models and eliminate the differences in the GEBVs associated with the form of prediction, the SG of each trait was estimated by the corrected phenotypic means and heritability from each model for the populations selected for the next selection cycle and seed donors. In this sense, the G-BLUP ST was chosen for the selection of juçaizeiro genotypes, because this model provided higher SG, in addition to being a simpler and more usual method for breeding programs.
The production of fruits per plant (PFP) is a feature that serves both the productive sector and the industrial sector, as it provides greater economic return to producers and increases the amount of raw material for industries that currently suffer from a shortage of this base product for their operations. In this sense, the selection differential amounts to an increase of 98.40%, which indicates a great advance for the culture due to the possibility of achieving great gains in productivity, resulting in an incentive to install of culture.
As the juçaizeiro is a wild and perennial species, which presents a series of characteristics intrinsic to its development, the practical activities of conventional breeding techniques become more difficult. In general, when evaluating the structure of a breeding program aimed at juçaizeiro, it is expected that the gains obtained per unit of time will be reduced. In this sense, the present study stands out for applying genomic prediction techniques in Euterpe edulis, a wild open-pollinated species, aiming to increase gains per unit of time.
Therefore, perhaps the biggest advantage brought by this methodology for a wild breeding population, as is the case of Euterpe edulis, is the use of the genomic matrix to bring knowledge of the genetic relationship between individuals for statistical analysis, which was previously unknown. Thus, the establishment of this information provides more accurate estimates and, consequently, greater reliability to the program. Furthermore, the use of genomic prediction may allow a reduction in the selection cycle by being able to predict the behavior of non-evaluated genotypes.
All the traits evaluated in the study influence the economic potential of the species, covering the interests of the productive sector and the industrial sector. Additionally, in function of the alterations observed by the selection of the six superior genotypes, we can conclude that the genomic prediction using the G-BLUP ST was efficient to provide alterations of the means in the desired directions. That is, obtaining an increase in the average phenotypic response of the selected genotypes for the characters of NB, MFB, PFP and PY, and a reduction for EDF.

Conclusion
Our results showed that the GBLUP-ST genomic prediction was more efficient in selecting the best genotypes, the selection provided substantial gains in the desired direction for multiple traits. Thus, the six selected genotypes (UFES.A.RN.390, UFES.A.RN.386, UFES.A.RN.080, UFES.A.RN.383, UFES.S.RN.098 e UFES.S.RN.093) can be used as commercial seed donor genotypes for the development of seedlings for the implantation of productive orchards, which will meet the demands of the productive sector and the consumer market.